/*----------------------------------------------------
Class: plant_2DSOS_RDM
DESCRIPTION:
This class implements random deposition(RD) with migration(M). The lattice
is 2D solid-on-solid(SOS).

REFERENCE/BASED ON:
n_sos_2D_L200_T650.c by Gangshi and Jianqiao
sos_2d_L100_T700.c

INPUT:
* T: Temperature, K
* W: deposition rate, layer/s
OUTPUT:
* z1~z4:
* R2: roughness square, layer^2
* M2: RMS slope square, dimensionless
* H: average height, layer
STATE:
* lattice: a 2D array save the height profile.
OPTIONS:
*

AUTHOR: Xinyu Zhang (zxy1256@gmail.com)
----------------------------------------------------*/
#ifndef PLANT_2DSOS_RDM_H
#define PLANT_2DSOS_RDM_H

#include "Signal.h"

const double EV2KB = 11608.696;
const double KB2H  = 2.08144796e10;

class plant_2DSOS_RDM:public block
{
private:
  // Options
  int    LatticeSize;
  int    mode;
  double E_D;
  double E_S;
  double E_N;

  // Inputs
  double T;
  double W;

  // States
  int**  lattice;
  int**  lattice_stat;
  double *p_migr;
  double *p_desorp;
  double ads_ttl,migr_ttl,desorp_ttl,rate_ttl;
  double MIGR_RT[5];
  double DESORP_RT[5];
  double k_m2d;

  // Outputs
  double *cov_z1;
  double *cov_z2;
  double *cov_z3;
  double *cov_z4;
  double h;
  double r2;
  double m2;

  // Others
  double **ef_cos_table; // A precalculated table for the cos part in eigenfunctions
  double **ef_sin_table; // A precalculated table for the sin part in eigenfunctions

  void add_remove_atom(int x, int y, int num);
  void adsorption();
  void migration();
  void desorption();
  double OneMCTime(double r_d);
  void migr_site(int *migr_x,int *migr_y);
  void desorp_site(int *desorp_x,int *migr_y);
  void set_ads_rate(double W_set);
  void set_temperature(double T_set);
  void GetNeighbors(int x,int y, int *xpos,int *ypos);	
  void init_migr_rate();
  void init_desorp_rate();
  void init_p_migr();
  void init_p_desorp();
  void init_local_stat(int x, int y);
  void h2mode();
  void h2mode_fftw();

  inline int index2D_h(int i,int j){return i*LatticeSize+j;};
  inline int index2D_mode(int m,int n){return m*(mode+1)+n;};

  // Eigen-functions
  inline double phi1mn(int m,int n,int ix,int iy){
    //return 2.0*sin(2*m*pi/LatticeSize*ix)*sin(2*n*pi/LatticeSize*iy)/pi;
    return ef_sin_table[m][ix]*ef_cos_table[n][iy];
  };

  inline double phi2mn(int m,int n,int ix,int iy){
    /*double coeff;
    if((m == 0 || m*2==LatticeSize)&&(n==0 || n*2==LatticeSize)){
    coeff = 1.0/pi;
    }
    else if(m!=0 && 2*m!=LatticeSize && n!=0 && 2*n!=LatticeSize){
    coeff = 2.0/pi;
    }
    else{
    coeff = SQRT_TWO/pi;
    }
    return coeff*cos(2*m*pi/LatticeSize*ix)*cos(2*n*pi/LatticeSize*iy);*/
    return ef_cos_table[m][ix]*ef_sin_table[n][iy];
  };

  inline double phi3mn(int m,int n,int ix,int iy){
    /*double coeff;
    if(n==0||n*2 == LatticeSize){
    coeff = SQRT_TWO/pi;
    }else{
    coeff = 2.0/pi;
    }
    return coeff*sin(2*m*pi/LatticeSize*ix)*cos(2*n*pi/LatticeSize*iy);*/
    return ef_sin_table[m][ix]*ef_sin_table[n][iy];
  };

  inline double phi4mn(int m,int n,int ix,int iy){
    /*double coeff;
    if(m==0||2*m==LatticeSize){
    coeff = SQRT_TWO/pi;
    }else{
    coeff = 2.0/pi;
    }
    return coeff*cos(2*m*pi/LatticeSize*ix)*sin(2*n*pi/LatticeSize*iy);*/
    return ef_cos_table[m][ix]*ef_sin_table[n][iy];
  };

  // The following declaration is to prevent compiler default implementation.
  plant_2DSOS_RDM(const plant_2DSOS_RDM&);
  plant_2DSOS_RDM& operator=(const plant_2DSOS_RDM&);
public:
  plant_2DSOS_RDM(double idt,int iLx=100,int iMode = 50,double iE_D = 1.8,double iE_S = 1.2,double iE_N=0.6);
  plant_2DSOS_RDM(const string FileName);
  ~plant_2DSOS_RDM();
  void update(double sysTime,double iT,double iW);
  virtual void update(double sysTime);
  virtual void reset();
  void   getSnapshot(const string FileName);
  double getM2();
  double getR2();
  double getH();
  double getTime();
};

double ranO(int *idum);
double RandRand();
double RandN1();
#endif
